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A rigorous proof for convergence of the Wolf method ^] for calculating electrostatic energy of a 
periodic lattice is presented. In particular, we show that for an arbitrary lattice of unit cells, the 
lattice sum obtained via Wolf method converges to the one obtained via Ewald method. 
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I. INTRODUCTION 
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The classical Madelung problem yj| has an important 
role in atomic and molecular simulations involving elec- 
trostatic interactions. Consider an arbitrary lattice with 
a unit cell that is composed of N charges {qi,...,qN} 
and let linearly independent vectors 81,82,63 G M"^ de- 
note the lattice vectors. We assume the charge neutrality 
condition for the unit cell, i.e., X]t=i 9* ~ ^- Then the 
Madelung problem for calculation of the total electro- 
static energy of the unit cell located at the origin can be 
expressed as 
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(the matrix with lattice 
denotes 



- Vj , where r,; 



where V = [ei 82 83] G 
vectors as its columns) and r^j — Yi 
the atomic position within the unit cell. The prime on the 
summation emphasizes that we exclude self-energy, i.e., 
for n = the term i = j is omitted. In order to be able to 
use the well-established theory of multi-dimensional zeta 
functions [2], we introduce the following non-standard 
representation of electrostatic energy of a unit cell in an 
arbitrary lattice: 
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where T stands for matrix transport and Qij = 2e,;-ej is a 
positive-definite matrix (twice the metric tensor), = 
Urjj with U = and the prime on the summation 

denotes the exclusion of the self energy. For example, for 
an orthorhombic lattice with lattice parameters a = |8i|, 
b ~ I82I and c = [83 1 we have 
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where a and c are unit cell parameters, we have 
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One method of calculating the above lattice sum is to 
use direct sums. However, it is a well known fact that 
([2]) is a conditionally convergent series, which means that 
^ is meaningless unless the order of summation of the 
terms is specified. It is interesting to note that sum- 
mation over regions that may seem to be natural can 
diverge. As an example, for a 3-dimensional NaCl-type 
ionic crystal, it was shown that this lattice sum does 
not converge over expanding spheres Q, expanding el- 
lipsoids, and some specific expanding polygons but 
it would converge for expanding cubes [3|. The direct 
summation method is not practical due to the slow rate 
of convergence. 

Another method for dealing with ^ is to find some 
analytic continuation for this expression over the com- 
plex plane and then to find some fast converging series 
to evaluate this analytic continuation. The celebrated 
Ewald method [Hi uses this procedure. One can write 
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with 



^q(s,p)= E' 
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where again the prime denotes the exclusion of any infi- 
nite summands. Zq{s, p) is a special case of the general 
Epstein zeta functions [2, It is known that ^q(s, p) 
is uniformly and absolutely convergent for any complex 
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number s with ?le(s) > 3 and it has a meromorphic con- 
tinuation to the whole complex s-plane 0, H, 01 ■ From 
here on, by Zq(s, p) we mean its analytic continuation. 
Zq(s, p) is analytic everywhere except for the simple pole 
at s = 3 with the residue 
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Using some special functions, it is possible to write 
Zq{s,p) and thus Sceii in terms of a rapid converging 
series that recovers the standard Ewald method [1, Q . 

Determining the relation between the above-mentioned 
two methods is not a straightforward task. Specifically, 
in what order one should sum the terms of the series 
([2]) to obtain the Ewald's result? Borwein, et al. Q 
showed that for NaCl-type crystals, summing over cubes 
would yield the Ewald result. It is interesting to note 



all of the summands are positive, partial sums of this 
series would be unbounded. Thus, it is meaningless to 
expect to find Zq(1,p) using direct sums. 

As we mentioned earlier, for NaCl-type crystals direct 
summation over expanding cubes converges while sum- 
mation over expanding spheres diverges. One may guess 
that what makes the expanding cubes to converge is that, 
unlike expanding spheres, each cube is charge neutral, 
and thus it may be possible to obtain a converging se- 
quence over spheres if one somehow converts the regular 
spheres to charge neutral ones. This is the main idea 
of the Wolf method @ and the earlier work of Buhler 
and Crandall In particular, for a general lattice of 
charges. Wolf, et al. @ suggested that putting a mir- 
ror charge on the surface of sphere for each charge inside 
the sphere and neglecting the charges outside the sphere 
results in a convergent sequence that converges to the re- 
sult obtained via Ewald method. Although they verified 
their method by considering several numerical examples, 
they did not present a rigorous proof. 

In this paper, we present the missing proof of the con- 
vergence of Wolf's method for calculating electrostatic 
energy of an arbitrary lattice of charges. We should men- 
tion that Buhler and Crandall presented a proof for 
NaCl lattice with unit charges. Here we generalize their 
proof to arbitrary lattices. Note that Wolf, et al. [9] 
presented two differrent methods, namely, damped and 
undamped methods. Here by Wolf method we mean the 
undamped method. In §2 we present the proof and in §3 
we mention some concluding remarks and future direc- 
tions. 



II. PROOF OF THE CONVERGENCE OF 
WOLF'S METHOD 

We use the contour integral method of [10] to prove the 
convergence of Wolf's method. Let us first review the re- 
quired preliminaries. We consider the following analytic 



continuation of ([2]) 
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Note that i5(Q, 1) = ^Ewaid = ^ccii- Let ^'(a;) denote the 
function 
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Then, for R G R"*" , where R"*" is the set of the positive real 
numbers, define the spherically truncated (finite) sum 
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Note that £^fl(Q, s) roughly denotes the direct sum of ([2]) 
over spheres centered at charges in the unit cell located 
at the origin with radius R. From (|10l) one can see that 
if 



R = 7^(n,p^J) := 
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for some values of i, j and n, then the charges located 
on the surface of the sphere would have the additional 
weight of 1/2. But this would be irrelevant in our analysis 
because as we will see in the sequel, we are interested in 
the behavior of a continues function of i? as i? — oo. 
This can be expressed in term of i?fl(Q, s) if R does not 
take the discrete values 7?.(n,py) for wGl? and i,j = 
1, . . . , iV, and so we can exclude those values. 

In the Wolf method one considers spheres with radii 
R centered at each charge of the unit cell at the origin 
and puts a mirror charge on the surface of the sphere for 
each charge inside the sphere including center charges. 
Then energy of the unit cell, fwoifj is calculated using 
only charges located on and inside these spheres. This 
can be written as 



f Wolf = lim 
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But for a lattice with N charges in its unit cell, the term 
XiLi 9? is bounded, and thus we have 
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FIG. 1: The contour Qc.t in the complex s-plane. 



For c G Perron's formula states that [11 1 
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Finally, for any c,T E IR+ with c > 3 let a^^T be 
the contour in the complex s-plane depicted in Fig. [T] 
Also let St = {seC\0< fHe(s) < 3 and \3m{s)\ < T}. 
Then, the following lemma holds. 

Lemma. Let e E M+ . Then there exist c,T E R+ with 
c > 3, such that for any w E St the following relation 
holds as i? —)• oo; 
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Proof: Let s = a + it and la = where 
'■'^-"^ Zci{s,p)R^ 
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We are going to show that all of the above integrals are 
bounded. For Ii and I5 wc have |Zq(s,p)| = |Zq(c, p)|, 
IR^'l = R" and [s{s - w)\-^ = 0{t-^). Thus, these inte- 
grals are OiR^T'^). 

Using the standard methods of analytic number the- 
ory and the Phragmen-Lindelof theorem [ll|, [l^] , we find 
that \Zq{s,p)\ = 0(i=-3/2) for 3 - c < cr < c, and 



hence h and h are 0(i?=r=-^/2). If 3 < c < 7/2 and 
we can conclude that Ii, I2, h and I5 are 
0(i?') for any e > 0. 
Next, note that 
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= J2 |An-dr^=:CA(s,d), (18) 
where A is a matrix with nonzero determinant. Define 

AA(s,d)=Vd^7r-^/2 p|^0^^(^^j) (^g) 

Then, we have the following functional equation Q 

AA(s,d) = e-l'i|'*AB(3-s,0), (20) 

where B = A^^ with — T denoting the inverse transpose. 
Substituting ^ into ^ yields 



Ca(s, d) = , , , ^ /,V Cb(3 - .s, 0). 
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Now we replace s by 3 — s in /a and use ([T8|) and ((2T|) to 
write the integral over the upper half segment as 
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where in the last step, since c > 3, we use the fact that 
Cb(s, 0) is uniformly and absolutely convergent over the 
integration path. Stirling's formula states that for any 
fixed strip a<(T<P,ast^oo [ij] 

log [r(cr -I- it)] = 

(^a + it- ^ log(M) - M + ^ log(27r) + 0{t~^). 

(23) 

Using (p3l) . we conclude that to bound (|22|) one needs to 
bound 
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where to > is an arbitrary constant. But using the 
method of stationary phase ^3|], it can be shown that 
this integral is 0(logi?) [13] and so I3 is 0(i?') for e > 0. 
To summarize, we have proved that Ii, . . . , are 0{R'^) 
for e > and hence (fT6|) holds. □ 

Now we state the main result of this paper. 

Theorem. There exists T G M+ such that if s £ St 
then 
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In particular, setting s — 1 yields 
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Proof: Choose c,T G R+ in accordance with the pre- 
vious lemma and let w € St- Consider the integral 
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Note that since c > 3, the contour of the integral is in 
the region of the uniform and absolute convergence of 
Also c — 9{e{w) > 0, and hence using ([T5|) one can 
write (P7|) as 
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Then, we use (|28|) to write 
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or equivalently 
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Next we choose the contour ac,T as in FigH] Let /3 denote 
the rectangle with vertices 81,82, 83 and 84 (see Fig[T]). 
We have 
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One can use Cauchy's Residue Theorem to calculate the 
right side of ([5^ . Refereing to (f^ . it is evident that 
FwiQ, s) has simple poles at s = and s = w and may 
have a simple pole at s = 3. Since E{Q, s) is analytic at 
s = and s — w, we obtain 



ResF^(Q,s) = -w-^E{Q,0), 

s=0 
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To calculate the residue at s = 3 we use ((29)) and ([8]) to 
write 
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where we used the charge neutrality condition for the 
unit cell in the last step. Thus, we have shown that 
although Zq(s,p) has a simple pole at s = 3, charge 
neutrality implies that residue of Fuj^Q, s) at s = 3 van- 
ishes. Therefore, using Cauchy's Residue Theorem and 
equations ([M]). (IMl) . and (I55t we conclude that 



7^ I F^{q,s)ds ^ -w-^E{q,0) + w-^R"''E{Q,w). 



(36) 
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Substituting dST]) and §6^ into ([32]) results in 
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On the other hand, with the aid of (O, (HH), and (|29|) 
the left side of (|37p can be written as 
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for an arbitrary e > 0. Thus, ([57)1 becomes 
^i?XQ,"')-;^^fl(Q,0) = 
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£;(Q,0) + £;(Q,w)+i?-'"O(i?^). (39) 



Since £He(u') > 0, upon taking the limit of ([5^ as i? oo 
and replacing w with s, we obtain (j25p . This completes 
the proof. □ 



III. CONCLUDING REMARKS 

A rigorous proof for the Wolf method is given in this 
paper. However, there are still some open questions. As 
we mentioned earlier, we prove that the undamped Wolf 
method converges to the Ewald sum. But usually the 
undamped method converges very slowly and this makes 
it unfavorable in practice. To resolve this issue. Wolf, et 
al. Q modified their method and introduced the damped 
method. The electrostatic energy computed via damped 
method, S^^if' is 
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erfc(a;) = 1 — erf(x), erf(x) 



e"* dt. (42) 



Since lim erfc(a;) = 0, using (^5)) we obtain 
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Although damped method converges fast, it converges 
to values that depend on the damping parameter and 
crystal structure [a, ■ This means that the right side 
of (|43| converges to values that depend on a. Thus, it 
does not seem that one can prove a theorem counterpart 
to one presented here for the undamped method. But 
one may find an interval(s) for the values of the damping 
parameter to control the error of the damped method. 
Also it is interesting to note that the variation of the 
calculated energy versus the radius of the charge natural 
sphere depends on the crystal structure: it may have 
an oscillatory behavior as for NaCl cryst al [91 or non- 
oscillatory behavior as for PbTiOa crystal [Ijl . 

We saw that the charge neutralization idea works for 
spherical expanding domains. One may want to see if 
this idea would work for other (convex) domains as well. 
Answer to this question can help to justify the use of 
the Wolf method for other geometries like free surfaces, 
slabs and regions near crystal defects. Besides calculating 
energy. Wolf, et al. proposed a method for obtaining 
forces exerted on charges again without a rigorous proof. 
As force is a vector quantity, there is an ambiguity on 
how one should project mirror charges on the surface of 
the sphere. Since energy only depends on the distance 
between charges such ambiguity does not arise in the 
calculation of energy. Providing rigorous proofs for the 
above questions will be the subject of future work. 
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